This site uses cookies. By continuing to use this site you agree to our use of cookies. To find out more, see our Privacy and Cookies policy.

THE BIERMANN CATASTROPHE IN NUMERICAL MAGNETOHYDRODYNAMICS

, , , , , , and

Published 2015 March 19 © 2015. The American Astronomical Society. All rights reserved.
, , Citation Carlo Graziani et al 2015 ApJ 802 43 DOI 10.1088/0004-637X/802/1/43

0004-637X/802/1/43

ABSTRACT

The Biermann battery effect is frequently invoked in cosmic magnetogenesis and studied in high-energy density laboratory physics experiments. Generation of magnetic fields by the Biermann effect due to misaligned density and temperature gradients in smooth flow behind shocks is well known. We show that a Biermann-effect magnetic field is also generated within shocks. Direct implementation of the Biermann effect in MHD codes does not capture this physical process, and worse, it produces unphysical magnetic fields at shocks whose value does not converge with resolution. We show that this convergence breakdown is due to naive discretization, which fails to account for the fact that discretized irrotational vector fields have spurious solenoidal components that grow without bound near a discontinuity. We show that careful consideration of the kinetics of ion viscous shocks leads to a formulation of the Biermann effect that gives rise to a convergent algorithm. We note two novel physical effects: a resistive magnetic precursor, in which a Biermann-generated field in the shock "leaks" resistively upstream, and a thermal magnetic precursor, in which a field is generated by the Biermann effect ahead of the shock front owing to gradients created by the shock's electron thermal conduction precursor. Both effects appear to be potentially observable in experiments at laser facilities. We reexamine published studies of magnetogenesis in galaxy cluster formation and conclude that the simulations in question had inadequate resolution to reliably estimate the field generation rate. Corrected estimates suggest primordial field values in the range $B\sim {{10}^{-22}}$–10−19 G by z = 3.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

Dynamo theories of proto- and extragalactic primordial magnetic fields, which endeavor to explain how those fields achieved their current strength and structure, work by amplifying small initial seed fields by means of turbulent plasma motions (Kronberg 1994; Kulsrud & Zweibel 2008). However, the induction equation of resistive MHD,

Equation (1)

always admits the solution ${\boldsymbol{B}} ({\boldsymbol{x}} ,t)=0$. This simple observation poses a problem for the generation of the required seed fields, as they cannot be created in ideal MHD starting from a field-free state.

There have been several proposals for generating the required seed fields from mechanisms such as primordial phase transitions, or from processes occurring during inflation (see reviews in Widrow 2002; Widrow et al. 2012). The Biermann battery effect (Biermann 1950) provides another popular resolution of this problem (Kulsrud et al. 1997; Widrow 2002; Kulsrud & Zweibel 2008). The effect, which arises in consequence of the large difference in the electron and ion mass, is attributable to small-scale charge separation in the plasma. Pressure forces produce much larger accelerations of electrons than of ions, and the relative acceleration of the two components results in charge separation that must be balanced by an electric field

Equation (2)

where ne and Pe are the electron number density and pressure, respectively. Since this field is not, in general, irrotational, it can act as a source of magnetic field in the induction equation,

Equation (3)

generating a nonzero ${\boldsymbol{B}} $ from an initially unmagnetized state.

The Biermann battery effect has been successfully invoked in numerical simulations exploring the generation of seed fields in cosmological ionization fronts (Subramanian et al. 1994; Gnedin et al. 2000), protogalaxies (Davies & Widrow 2000), and Population III star formation (Xu et al. 2008). Magnetic field generation by the Biermann mechanism is also of significant interest in direct-drive and indirect-drive inertial confinement fusion, where strong gradients behind the converging shock can lead to dynamically important field strengths (Srinivasan & Tang 2013), and more generally in the field of high-energy density physics (HEDP), where the effects of field generation (Gregori et al. 2012) and amplification (Meinecke et al. 2014) can be examined in a laboratory setting at laser facilities, in experiments where large gradients are produced in strong plasma shocks (Fryxell et al. 2010; Tzeferacos et al. 2012, 2014).

While the generation of magnetic fields by the Biermann effect as a result of strong misaligned density and temperature gradients in the smooth flow behind shocks is well known, we show that there exists a previously unrecognized Biermann effect due to the electron–ion charge separation that occurs within ion viscous shocks. This effect is a consequence of the kinetic theory of shock structure in plasmas, as we show below.

A straightforward implementation of the Biermann effect in finite-volume Eulerian, purely Langrangian, and ALE codes, whether as a split source term or as a flux term, does not capture this physical process, and worse, it leads to nonconvergent results (Fatenejad et al. 2013). In symmetric situations such as planar or spherical shocks, where no field should arise, such implementations produce anomalous field generation that grows without bound with resolution (Fatenejad et al. 2013). This behavior is observed across a range of different MHD codes (see the discussion of codes in Fatenejad et al. 2013). This is the Biermann catastrophe of numerical MHD. We show that this failure is intimately related to the failure of such codes to correctly model the structure of the plasma shock.

In a gasdynamic/MHD formulation, where the shocks are modeled as zero-width discontinuities of the flow, the trouble arises from the behavior of the Biermann flux, which is to say from the electric field, Equation (2), in the vicinity of a shock. The gradient $\nabla {{P}_{e}}$, which, analytically speaking, acquires a Dirac δ component at the shock, is ascribed a numerical magnitude that grows without bound at the shock with increasing resolution. It is this divergence that is connected with failure of MHD codes to correctly predict shock-driven magnetic field generation in supposedly simple test cases, where the correct value of the generated field is zero. The investigation of this failure is one of the central concerns of this article.

That this failure was not previously recognized is a consequence of the history of the numerical study of the Biermann battery effect as a source of cosmic seed fields, which has a curious feature with respect to shocks. On the one hand, the importance of shocks to lifting barotropic constraints, thus making available the kind of nonaligned gradients required to drive the Biermann effect, has been widely recognized (Kulsrud et al. 1997; Davies & Widrow 2000). On the other hand, the direct effect of shocks—as opposed to that of their trailing downstream gradients—on this magnetogenesis has not really been carefully examined. This is probably because the computational strategies that have been adopted both circumvent difficulties at shocks and direct attention away from the magnetizing properties of shocks. For example, Kulsrud et al. (1997) and Gnedin et al. (2000) perform essentially hydrodynamic simulations, in which magnetic fields evolved by the induction equation have no dynamical role. Davies & Widrow (2000) forgo the induction equation altogether, in favor of the equation of inviscous evolution of vorticity ${\boldsymbol{ \omega }} $

Equation (4)

which, by comparison with Equation (3), implies that inviscous, nonresistive plasmas satisfy the equation

Equation (5)

(Kulsrud & Zweibel 2008). On its face, this relation would appear to suggest that in an initially unmagnetized and irrotational plasma, the magnetic field is simply proportional to the vorticity, ${\boldsymbol{B}} =\frac{c{{M}_{i}}}{e(1+\bar{Z})}{\boldsymbol{ \omega }} $, where Mi is the ion mass and $\bar{Z}$ is the average ionization fraction.

One difficulty with these approaches is that they entirely neglect to treat the field generation within the shock itself, as we noted above. There are indeed large, nonaligned gradients at the simulated shocks, but these gradients are unphysical side effects of the numerical strategies used to integrate the hydrodynamic equations, and consequently they are resolution dependent. It is therefore hopeless to expect convergence with resolution of the resulting magnetic fields.

Furthermore, any magnetization generated by the shock itself imprints itself on the magnetic field structure as the shock moves through, leaving behind a substantial residue that is superposed on the smooth-flow Biermann-generated field. It is essential to come to a correct understanding of the behavior of the Biermann term at shocks, to have any confidence that results arrived at in this manner bear any resemblance to reality. In particular, the presumption that shock jump condition on vorticity should bear a relation to the jump condition on magnetization that preserves the proportionality between the two has been hypothesized but never demonstrated (Kulsrud et al. 1997; Davies & Widrow 2000), and it seems in fact far-fetched: it would imply that magnetization is vorticity, that is, that the coincidence expressed by Equation (5) in fact reflects a deep identity, and that magnetic degrees of freedom of ideal MHD are somehow already contained in unmagnetized Eulerian hydrodynamics, encoded in derivatives of the velocity field. Such a claim is difficult to accept, and to even attempt to support it would require an analysis of the modification brought by the Biermann battery term to the Rankine–Hugoniot(R–H) jump condition on ${\boldsymbol{B}} ,$ an analysis that we furnish for the first time in this paper. We will return to a discussion of vorticity in Section 7.

There have also been full MHD simulations of magnetogenesis through the Biermann effect in galaxy clusters (Xu et al. 2008), which were performed using the ENZO code (Xu et al. 2010, 2011; Bryan et al. 2014). No convergence study of Biermann-generated magnetic fields was reported for these simulations, possibly because of the paucity of nontrivial analytical verification tests against which the code's implementation of the Biermann effect could be tested.

It is clearly long past time that the behavior of the Biermann effect at an MHD shock be fully analyzed and characterized, so as to furnish a mathematical-physics target at which algorithmic implementations can shoot. In order to do so, it is necessary to draw connections between the well-understood kinetic theory of planar plasma shocks (Shafranov 1957; Jaffrin & Probstein 1964; Zel'dovich & Raizer 2002) and that of spatially inhomogeneous MHD shocks, connections that have not so far been made or exploited in the literature.

In particular, it is essential to address the following questions.

  • 1.  
    How does the kinetic theory of the structure of plasma shocks inform our understanding of the Biermann effect in shocks, and how can it guide us to an accurate and convergent treatment of the effect in MHD codes?
  • 2.  
    Do we have any right to expect convergence with resolution from an MHD simulation with a source term as ostensibly misbehaved at shocks as that of Equation (2)? In other words, is the Biermann source term mathematically well defined near a discontinuity of a weak solution of the MHD equations, and if so, how does it affect the R–H jump condition on ${\boldsymbol{B}} $?
  • 3.  
    Assuming that the Biermann source term, Equation (2), can be interpreted in a sensible manner in the vicinity of an inviscid shock, are its predictions consistent with the predictions of kinetic theory near a shock? Or does the Biermann term in MHD need to be flux limited, in the style of thermal and radiation source terms whose misbehavior must be limited on kinetic theory grounds when gradients grow too large?

We address question (1) in Section 2, abstracting the essential ingredients of plasma shock theory required to construct a valid MHD model of the Biermann effect due to shocks. We address the remaining two questions in Section 3, where we demonstrate that (2) in fact the Biermann source of Equation (2) is mathematically consistent and well behaved near weak-solution discontinuities, and (3) the Biermann source term in fact yields the correct electromotive force (EMF) across the shock, matching the prediction of plasma shock theory (Zel'dovich & Raizer 2002; Amendt et al. 2009; Jaffrin & Probstein 1964).

We clarify the origin of the Biermann catastrophe as a numerical effect attributable to the difficulty of discretizing the source term of Equation (2) in the vicinity of a shock. We show that the numerical anomaly can be eliminated by leveraging the continuity of the electron temperature Te across shocks—a benefit of the kinetic-theory connection. Reformulation of the Biermann source term in terms of Te allows the singularity to be isolated and the flux of the magnetic field due to the Biermann effect to be rewritten in a manifestly finite form suitable for translation into a convergent numerical algorithm.

The Biermann effect is due to electron–ion charge separation and is sensitive to departure from thermal equilibrium between electrons and ions. Such a departure is precisely what occurs at shocks, so that a correct treatment of the effect at shocks necessarily requires that the disequilibrium be modeled. For this, a two-temperature plasma model is mandatory. An interesting consequence of this observation is the fact that the Biermann effect is enhanced not only at shocks but also at contact discontinuities, at ionization fronts, and at material species boundaries, because the electron partial pressure is discontinuous at such surfaces, even though the total pressure is continuous there. This enhancement cannot be computed in an equilibrium treatment of the plasma.

An additional point of this paper is to point out two novel and interesting effects associated with the Biermann effect in the neighborhood of a shock. In Section 4.1 we show that a resistive magnetic precursor is generated in resistive MHD, wherein the magnetic field generated by the Biermann effect in the shock "leaks" resistively from the shock into a region of the upstream fluid whose physical extent is proportional to the resistivity. In Section 4.2 we show that a thermal magnetic precursor is generated ahead of the shock through the Biermann effect by plasma motions generated by the shock's electron thermal conduction precursor. Both effects are potentially observable in laboratory conditions at high-intensity laser facilities such as Vulcan, Omega, and the National Ignition Facility (NIF). We discuss the observability of these effects at laser facilities. We show that appropriately designed experiments at such facilities could currently observe the resistive magnetic precursor, providing a clean experimental validation test of the Biermann effect in plasma shocks. We also argue that the smaller thermal magnetic precursor might become observable in future experiments.

2. REVIEW OF KINETIC THEORY OF PLASMA SHOCKS

We begin by reviewing some essential results from the kinetic theory of shocks in plasmas. The basic theory of the fluid structure of planar shocks in plasmas was set out in Shafranov (1957), while the electromagnetic structure of such shocks was discussed in Jaffrin & Probstein (1964) and more recently in Amendt et al. (2009). An extremely lucid presentation of these results may be found in Chapter 7 of Zel'dovich & Raizer (2002).

There are three essential ingredients to be imported from the kinetic theory of plasma shocks in order to fashion a working MHD model of the Biermann effect: the loss of thermal equilibrium between electrons and ions at shocks; the adiabatic behavior of electrons, up to electron heat conduction; and the charge-separation-induced electric field across the shock front. We now review these in order, and we conclude by presenting the MHD model that forms the basis for the rest of this paper.

2.1. Ion–Electron Disequilibrium

As discussed in Drake (2006), a strong shock disturbs the thermal equilibrium between electrons and ions in a plasma. That equilibrium is maintained by electron–ion collisions and operates over timescales ${{\tau }_{ei}}$ that are long compared to the shock-crossing time of a parcel of fluid entering the ion viscous shock in consequence of the large ratio ${{m}_{i}}/{{m}_{e}}$ (Spitzer 1962). As a result, it is essential to describe the fluid in terms of an additional degree of freedom—the electron temperature, Te—with respect to the usual equilibrium MHD model.

Fortunately, it is unnecessary to model the fluid using new inertial degrees of freedom to describe the electron fluid. A single inertial component for the fluid as a whole gives an adequate description of the fluid structure near the shock and a completely satisfactory description of the fluid in smooth flow regions, where electron–ion collisions restore local thermal equilibrium (Drake 2006; Zel'dovich & Raizer 2002). This makes it easier to adapt existing MHD codes to treat the Biermann effect correctly, since it is much easier to add a scalar degree of freedom for Te than it would be to deal with two velocity fields, one each for the electron fluid and for the ion fluid.

2.2. Nearly Adiabatic Electrons

The large ion-to-electron mass ratio, which, as we recall from the discussion preceding Equation (2), is responsible for the charge separation that produces the Biermann effect, has the further consequence that electron mobility is much higher than ion mobility, and that thermal conductivity due to electron-electron collisions dominates the heat transport in the fluid. At the same time, the forces between electrons and ions as the fluid crosses the ion viscous shock front are effectively dissipation-free, as they are simply electrostatic fields generated by charge separation, and collisional dissipation processes are too slow to act during the shock-crossing timescale.

These observations lead to the conclusion that while the ions undergoing shock compression experience the usual thermodynamically irreversible, entropy-generating process associated with shocks, the electrons are compressed adiabatically by the electrostatic forces exerted upon them by the ions and consequently do not suffer entropy increments due to shock compression. The electron entropy would be a passively advected scalar, then, except for the dissipative effect of electron thermal conduction and for the slow (relative to timescales relevant to the shock) effect of electron–ion collisional equilibration. Electron thermal conduction effectively rules out any sudden change in Te at the shock, since such a change would produce an enormous restoring heat flux to heal the discontinuity.

The fluid structure that follows from these considerations is of a sudden discontinuous compression of the ions at the shock, accompanied by a smooth increase in Te, which is continuous throughout the shock (in the Eulerian limit where the shock width tends to zero, Te acquires a discontinuous derivative in the direction of the shock normal). The electron temperature exhibits a thermal precursor that leads the shock. The size ${{\lambda }_{T}}$ of the precursor region may be calculated by balancing advection against heat diffusion in the frame of the shock and is about

Equation (6)

where ${{\kappa }_{e}}$ is the electron thermal conductivity, ${{c}_{v,e}}$ is the electron specific heat at constant volume per unit mass, and D is the shock speed (Shafranov 1957; Zel'dovich & Raizer 2002). The effect of this precursor region on accretion shocks in galaxy clusters has been recently studied in Smith et al. (2013).

2.3. Electric Structure

The sudden change in the motion of ions entering the shock sheath, together with the more highly mobile motion of the electrons in the vicinity of the shock, results in a zone of charge-separation-driven electric field in the normal direction to the shock (Jaffrin & Probstein 1964; Amendt et al. 2009; Zel'dovich & Raizer 2002). By solving the Navier–Stokes equation across the shock (treating the results as means of kinetic-theory distributions, since a fluid picture is clearly invalid inside the viscous shock sheath), Jaffrin & Probstein (1964) calculated the full electric structure across the shock. We do not require the full structure here, since in the end we wish to adopt an Eulerian picture of the shock, in which the width of the shock is zero. For our purposes, the EMF across the shock is all that is required. It is shown in Amendt et al. (2009) that the EMF $\mathcal{E}$ is given by

Equation (7)

where the subscripts "u" and "d" denote "upstream" and "downstream," respectively, and where we have assumed complete ionization to obtain the second line.

This electric field is due to charge separation and thus has a common ancestry with the electric field ${{{\boldsymbol{E}} }_{B}}$ in Equation (2). However, it is proper to the shock, rather than to the smooth-flow portion of the fluid. It is essential, for the purpose of physical consistency, to demonstrate that an MHD model implementing the Biermann effect should reproduce the shock-crossing EMF given in Equation (7). That this is possible is demonstrated in the next section.

2.4. The MHD Model

We wrap up this section by exhibiting the Biermann-effect-laden MHD model implied by these kinetic-theory considerations. The dynamical system describing the model is given in conservation form as follows:

Equation (8)

Equation (9)

Equation (10)

Equation (11)

Equation (12)

In these equations, $P={{P}_{e}}+{{P}_{i}}$ is the total material pressure, ${{\epsilon }_{T}}$ is the total specific thermal energy per unit mass, and se is the specific electron entropy per unit mass.

Equations (8)–(9) are the usual MHD equations of mass and momentum conservation. The energy conservation equation, Equation (10), includes a term to account for resistive dissipation and another to account for the Biermann effect.

Equation (11) expresses the conservation of electron entropy se (up to heat conduction and electron–ion heat exchange), which, as we pointed out above, is required by the same approximation ${{m}_{e}}/{{m}_{i}}\to 0$ that gives rise to the Biermann effect in the first place. This equation is introduced to describe the additional degree of freedom required by the two-temperature plasma treatment.

Finally, in the induction equation, Equation (12), the Biermann effect is expressed here as a curl, rather than in the conservation form—a divergence of a flux tensor—that it will assume later.

These dynamical equations are supplemented here by perfect gas equations of state (EOSs) for both the electrons and the ions. We assume total ionization throughout what follows for the sake of simplicity.

3. THE BIERMANN EFFECT AT SHOCKS

3.1. Kinematics of the Biermann Effect at Shocks

By inspection of the Biermann source term in Equation (3), we see that it is proportional to $\nabla {{n}_{e}}\times \nabla {{P}_{e}}$. It follows that field generation by the Biermann effect can only occur if the gradients of ne and Pe are not aligned. This means that in shocks with planar, cylindrical, or spherical symmetry, the field generation rate is zero. We must therefore treat nonsymmetric shock surfaces in order to analyze nontrivial cases of field generation. Accordingly, in this subsection we set up the basic kinematics of the Biermann effect at general shock surfaces.

To describe the shock surface, we introduce a level function ${\Psi }({\boldsymbol{x}} ,t)$ and use it to define the shock surface ${\Psi }({\boldsymbol{x}} ,t)=0$. Note that the function Ψ is of no dynamical significance, but is rather simply a mathematical convenience for describing the shock surface. We will assume that the shock is moving in the direction of the normal vector ${\boldsymbol{n}} \equiv \nabla {\Psi }/|\nabla {\Psi }|$, so that the region ${\Psi }({\boldsymbol{x}} ,t)\gt 0$ is upstream, whereas the region ${\Psi }({\boldsymbol{x}} ,t)\lt 0$ is downstream. These definitions are illustrated in Figure 1.

Figure 1.

Figure 1. Illustration of the kinematics of the Biermann effect at a shock surface.

Standard image High-resolution image

We denote the local shock speed along ${\boldsymbol{n}} $ by D. By considering the motion of the level surface ${\Psi }({\boldsymbol{x}} ,t)=0,$ it is not difficult to show that

Equation (13)

To describe a moving MHD discontinuity that coincides with the surface ${\Psi }({\boldsymbol{x}} ,t)=0$, we will frequently express a field quantity X by the decomposition $X={{X}_{u}}({\boldsymbol{x}} ,t){\Theta }({\Psi })+{{X}_{d}}({\boldsymbol{x}} ,t){\Theta }(-{\Psi })$, where Xu and Xd are continuous functions, and where ${\Theta }({\Psi })$ is a Heaviside step function. After substituting such expressions into field equations, we will find some terms proportional to the Dirac distribution $\delta ({\Psi })$, resulting from differentiation of the Heaviside functions. We will refer to the collected terms of this form as the "shock part" of the evolution equations. As will be seen, such terms embody R–H jumps, which can thus be efficiently extracted for these nonsymmetric shocks. This trick is also useful for deducing hydrodynamic flux terms, as will also be evident below.

It is convenient to reformulate the Biermann source term using the ideal EOS to replace ne with Te, the electron temperature. Assuming an ideal gas EOS, the reformulated electric field is

Equation (14)

so that the source term due to the Biermann effect in the induction equation is

Equation (15)

An important reason for this reformulation is that, as we saw in Section 2, in the presence of electron thermal conduction, Te is continuous at the shock, whereas ne is not (Shafranov 1957; Zel'dovich & Raizer 2002). As we saw, the continuity of Te is a consequence of the high mobility of electrons relative to ions and is therefore part and parcel of the same approximation that led to the Biermann source term (Equation (2)) in the first place. It is central to the developments that follow.

The log electron pressure ${\rm ln} {{P}_{e}}$ is discontinuous at the shock and may be represented at time t by

Equation (16)

where ${{P}_{e,u}}({\boldsymbol{x}} ,t)$ and ${{P}_{e,d}}({\boldsymbol{x}} ,t)$ are continuous functions.

The discontinuity of Pe leads to a Dirac δ singularity in $\nabla {{P}_{e}}$. Using the distributional relation $d{\Theta }({\Psi })/d{\Psi }=\delta ({\Psi })$, we have

Equation (17)

We may easily interpret the term in Equation (17) that is proportional to $\delta ({\Psi })$ as the gradient proper to the shock, and the remaining terms as the gradient in smooth flow. Since we are interested in the shock behavior, we define

Equation (18)

where we have used the continuity of Te and the assumption of complete ionization to replace the pressure ratio across the shock with the corresponding compression ratio ${{\rho }_{u}}/{{\rho }_{d}}$.

The electron temperature Te is continuous across the shock, but its normal derivative is discontinuous. We may therefore represent Te by

Equation (19)

where ${{T}_{e,u}}({\boldsymbol{x}} ,t)$ and ${{T}_{e,d}}({\boldsymbol{x}} ,t)$ are continuous functions, and where continuity at the shock implies

Equation (20)

By virtue of the continuity of Te, the gradient $\nabla {{T}_{e}}$ is not burdened by a Dirac δ, and we obtain

Equation (21)

Inserting Equations (18) and (21) into Equation (15), we obtain

Equation (22)

where we have made use of the distributional relation

Equation (23)

where, to get the third line, we have used the fact that the Heaviside function squares to itself, since its values are either 0 or 1. This derivation of ${\Theta }\delta =\delta /2$ is admittedly cavalier in its treatment of distributional quantities and is presented for brevity. The result may also be obtained by a limiting procedure, in which the Heaviside function is represented by the limit of a family of continuous functions—this is, after all, the "Eulerian" limit of the description of a shock as the viscosity goes to zero.

Introducing the shock normal vector ${\boldsymbol{n}} $, and using the fact that the tangential derivative of Te is continuous (so that ${\boldsymbol{n}} \times \nabla {{T}_{e,u}}={\boldsymbol{n}} \times \nabla {{T}_{e,d}}$ at the shock), we obtain

Equation (24)

Equation (24) is amenable to a clarifying interpretation, which follows from the recognition that if V is a region containing some portion Σ of the shock surface, and the differential element of surface area on Σ is dA, then ${{\int }_{V}}{{d}^{3}}{\boldsymbol{x}} \;\delta ({\Psi }({\boldsymbol{x}} ))|\nabla {\Psi }|={{\int }_{{\Sigma }}}dA$, which is to say that $\delta ({\Psi }({\boldsymbol{x}} ))|\nabla {\Psi }|{{d}^{3}}{\boldsymbol{x}} $ is the differential area element on the shock surface. We may therefore interpret the quantity $(c{{k}_{B}}/e){\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}} $ in Equation (24) as a field generation rate per unit time per unit area of the shock surface. This interpretation already allows us to see that the Biermann source term should give rise to finite, mathematically sensible field generation.

To verify this, we may now obtain the field generation rate due to the passage of the shock by inserting the above expressions into the induction equation, Equation (3), neglecting resistivity. The equation may be written as

Equation (25)

In the neighborhood of the shock, we define

Equation (26)

Equation (27)

Equation (28)

Equation (29)

where the "tangential" components satisfy ${{{\boldsymbol{B}} }_{Tu}}\cdot {\boldsymbol{n}} =0$, ${{{\boldsymbol{B}} }_{Td}}\cdot {\boldsymbol{n}} =0$, ${{{\boldsymbol{u}} }_{Tu}}\cdot {\boldsymbol{n}} =0$, ${{{\boldsymbol{u}} }_{Td}}\cdot {\boldsymbol{n}} =0$. By virtue of the solenoidal condition $\nabla \cdot {\boldsymbol{B}} =0$, Bn is continuous across the shock.

We transform the advection term $\nabla \cdot \left( {\boldsymbol{Bu}} -{\boldsymbol{uB}} \right)$ in Equation (25) by inserting Equations (26)–(27), grouping the Heaviside functions, and collecting the "shock" terms proportional to the Dirac δ that arises from differentiating the Heaviside functions. Using Equation (13) to replace $\partial {\Psi }/\partial t$, and using Equation (24), we obtain

Equation (30)

Equation (31)

hence the R–H jump condition, generalized by the Biermann flux, is

Equation (32)

Here the notation $[\ldots ]_{d}^{u}$ means the difference of the upstream and downstream values at the shock location. The first two terms in Equation (32) constitute the usual jump condition for the induction equation (see, e.g., Chapter 7 of Gurnett & Bhattacharjee 2005). The final term is the contribution to the jump condition from the Biermann effect, which is seen to produce a finite, well-defined discontinuity at the shock. We may obtain a useful result by considering the jump in ${\boldsymbol{B}} $ due to a shock advancing into a quiescent, unmagnetized fluid (Bn = 0, ${{B}_{Tu}}=0$, ${{{\boldsymbol{u}} }_{u}}=0$). The downstream magnetic field is then given by

Equation (33)

where we have used the condition of mass continuity at the shock, ${{\rho }_{u}}D={{\rho }_{d}}(D-{{u}_{nd}})$, to replace the term $(D-{{u}_{nd}})$, and where we have also used the continuity of Te to replace the pressure ratio in the log with the density ratio.

We conclude from the above development that the Biermann source term is mathematically well defined even at weak solution discontinuities and yields definite finite predictions to which a properly designed numerical algorithm should be expected to converge.

3.2. Physical Compatibility of the Biermann Source Term with Plasma Shock Theory

In the Introduction, we raised the question of whether the Biermann source term, Equation (2), behaves near a shock according to the predictions of kinetic theory, as summarized in Section 2. We now explain the line of thinking behind this question and answer it definitively.

The induction equation, Equation (3), can be cast in conservative form. In this form, the Biermann source term assumes the form of the divergence of a flux tensor, whose components are linear in $\nabla {{P}_{e}}$. It is clear that the construction of these flux components from the gradient of a discontinuous function is in some way associated with the numerical troubles that arise with the Biermann effect near shocks.

The key point here is that "trouble with a flux computed from a derivative" is actually a familiar situation from radiation diffusion, where the radiation flux, which is proportional to the gradient of the radiation temperature, yields unphysical (superluminal) fluxes at regions where the temperature changes sharply (see, e.g., Mihalas & Weibel-Mihalas 1999). Another analogous situation arises with respect to electron thermal conduction, where the thermal flux from $\nabla {{T}_{e}}$ may yield unphysically large transport at shocks, in consequence of the discontinuous change in Te (Mihalas & Weibel-Mihalas 1999). In both of these cases, a straightforward work-around is furnished by a flux limiter, in which a maximum flux deduced from a kinetic-theory picture of the transport is used to smoothly cut off the flux in regions where gradients get large, while leaving the fluxes in smooth regions unaltered.

This is the reason that the question of the validity of the Biermann flux at shocks is a natural one to ask. If it were the case that the flux is just wrong when $|\nabla {{P}_{e}}|$ gets too large, a reasonable approach would be to use the limiting value for the flux from kinetic theory (Equation 7) and impose that flux in the shock as a cutoff through some kind of flux limiter. If the electric field due to charge separation given by Equation (2) results in an EMF across a shock that exceeds this value, we could cut it off at this value and obtain approximate results analogous to flux-limited approximation-based results from diffusion-limited radiation transport.

Being thus prepared for bad news about the Biermann term from kinetic theory, we instead are met by a surprise upon examination of the MHD flux. The Biermann electric field at a shock is, according to Equations (14) and (18),

Equation (34)

Integrating this over a vanishingly short shock-crossing path l along the normal to the shock, we obtain the electromagnetic work done by an electron in the presence of the Biermann electric field,

Equation (35)

which is the same as the value inferred from kinetic theory, Equation (7).

We conclude from this argument that there is no analogy to flux-limited diffusion in the behavior of the Biermann source term at shocks. The "plain vanilla" source gives the correct EMF across the shock and needs no flux limiter to cut it off there. It should be perfectly possible to construct a physically valid algorithm to represent the Biermann MHD source term unmoderated by limiters, one that faithfully reproduces the predictions of kinetic theory at a shock.

3.3. Origin of the Numerical Biermann Catastrophe

As discussed in the Introduction, MHD simulations incorporating the Biermann battery effect have invariably produced catastrophic nonconvergent behavior as soon as discontinuities develop. There are two (related) kinds of catastrophes that have manifested themselves: simulations with spherical or planar symmetry, which develop a nonzero field at shocks with a field intensity that grows with increasing resolution, despite the fact that the charge-separation electric field is irrotational and the magnetic field generation rate should therefore be zero (Fatenejad et al. 2013); and simulations of HEDP laser experiments containing spatially inhomogeneous shocks, similar to the (highly simplified) simulations presented below in Section 6.2.2, in which the field generated at shocks, which is not expected to be zero, fails to converge to a finite value with increasing resolution—in order to obtain reasonable results at all in such simulations, the cells participating in shocks, contact, and material discontinuities must be located at each time step, and the generation rate in those cells artificially set to zero (Fatenejad et al. 2013). Such simulation results are clearly not correct, since they do not treat magnetogeneration due to the ion viscous shock, but at least they do converge with resolution.

Clearly something goes very wrong with the numerics when the flow develops discontinuous behavior, whereas the behavior of the Biermann term in smooth flow is finite, convergent, and stable. To investigate the origin of these catastrophes, we have addressed the following question: suppose the electric field ${{{\boldsymbol{E}} }_{B}}$ of Equation (15) is irrotational, in the sense that

Equation (36)

To what extent is this irrotational character preserved after Te and Pe are separately discretized in a volume-based scheme?

We consider the discretization of two scalar functions $f({\boldsymbol{x}} )$ and $g({\boldsymbol{x}} )$ that are assumed to have gradients that are everywhere collinear, so that there exists some function $\alpha ({\boldsymbol{x}} )$ such that $\nabla g=\alpha \nabla f$. A pair of such functions satisfy $\nabla g\times \nabla f=0$, so that the vector field $g\nabla f$ is irrotational. The Taylor series of such a pair of functions are related to each other by the collinearity constraint. We Taylor-expand both functions to third order, average the resulting expansion over a grid of control volumes, and compute the "discrete curl"—that is, the circulation integral about a cell face—of the discretized "$g\nabla f$" vector field. We omit the details of this extremely lengthy calculation for the sake of brevity and merely present the conclusions that may be drawn therefrom.

The result of the calculation is that the leading-order behavior of the discrete curl of such a discretized vector field is $\mathcal{O}({{{\Delta }}^{2}})$, where Δ is the grid spacing. The coefficient multiplying this term is homogeneous of the order of 3 in the derivatives of f and g. What this means is that in regions of space where the two functions are smooth (C2 or smoother), the coefficient has a finite approximation in the discretization, and the discrete curl tends to zero as ${{{\Delta }}^{2}}$ with increasingly fine resolution, yielding the correct curl of an irrotational vector field in the limit.

On the other hand, in the presence of a discontinuity, the coefficient of the discrete curl has no finite approximation in the discretization, but rather diverges as ${{{\Delta }}^{-3}}$ when ${\Delta }\to 0$, reflecting the meaninglessness of the Taylor approximation in proximity to a discontinuity. As a consequence, the discrete curl as a whole diverges as ${{{\Delta }}^{-1}}$ with vanishing Δ in the neighborhood of a discontinuity.

Setting $f={\rm ln} {{P}_{e}}$ and $g={{T}_{e}}$, we recognize immediately the source of the discrete pathology described above. The discretization of Pe and Te produces a spurious solenoidal component of ${{{\boldsymbol{E}} }_{B}}$ above and beyond any real, physically correct solenoidal component. This solenoidal anomaly is small and converges to zero with increasing resolution in smooth flow. Near a discontinuity, however, the anomaly is not convergent, but rather diverges as ${{{\Delta }}^{-1}}$.

This is the explanation of the numerical Biermann catastrophe. It is, fundamentally, a completely predictable failure of naive, stencil-based approximations to the Biermann source term (Equation (2)), which are not meaningful when Taylor-series approximations to the fluid variables fail in the presence of a discontinuity. The failure is irreducible and does not depend on whether the Biermann flux is added to other MHD fluxes or computed separately as a split term.

3.4. The Cure for the Biermann Catastrophe

The above diagnosis of the origin of the numerical Biermann catastrophe does not directly suggest a cure for the problem. That a cure should exist, however, is strongly suggested by the fact that the analytic treatment of the Biermann effect at shocks in Section 3.1 leads to finite and physically sensible field generation rates. We therefore now partly retrace that development with a view to casting those results in a form suitable for formulation as a numerical algorithm to be incorporated in a Godunov-type volume-based MHD scheme.

In order to clarify this proposed strategy, let us briefly review how volume-based hydrodynamics/MHD schemes work (for details, see, e.g., LeVeque 2002). Nonlinear systems of hyperbolic partial differential equations (PDEs) such as MHD are prone to develop discontinuities such as shocks and contacts. Direct finite-difference discretization approaches break down when this occurs. This breakdown is circumvented by appealing to the conservative nature of the equations and replacing the differential equations with their corresponding conservation-integral forms, applied to a regular mesh of control volumes. Field quantities are interpreted as volume averages over control volumes and are updated by time-centered (i.e., half-time-step forward) fluxes across control-volume faces. These fluxes are finite even when solution discontinuities cause the corresponding differential PDE source terms to misbehave.

The crucial procedure for carrying out such schemes successfully is the computation of the time-centered fluxes by means of the solution of Riemann shock tube problems at the cell interfaces. A Riemann solver takes as its input piecewise-constant initial data with a discontinuity at the interface and calculates the resulting propagating state, comprising piecewise-constant regions separated by discontinuity waves and rarefactions belonging to different families of characteristics. This propagating state is used to infer the MHD variables at the interface a half time step following the initial state. MHD fluxes computed from these variables are then used to update the cell-averaged field quantities (LeVeque 2002; Toro 2009).

Riemann solvers relate adjoining MHD states separated by a discontinuity traveling with wave speed D by means of the R–H relations,

Equation (37)

Here ${\bf \Phi }$ is the state vector of conserved field quantities, ${{{\bf \Phi }}^{T}}\equiv [\rho ,\rho {{{\boldsymbol{u}} }^{T}},\rho \mathcal{E},{{{\boldsymbol{B}} }^{T}}]$, that is, the mass, momentum, and energy densities and the magnetic field. ${\boldsymbol{F}} $ denotes the vector of fluxes corresponding to ${\bf \Phi }$. The subscripts "L" and "R" denote "left" and "right" states, respectively, and D is positive for a wave traveling from left to right. (Toro 2009).

We can now state a minimal requirement in order to assimilate the Biermann flux to the other MHD fluxes used to update the fluid state: We need to establish the terms added by the Biermann effect to the flux vector ${\boldsymbol{F}} $. Once these additive fluxes are determined, they can be added to the ideal MHD fluxes.

The additional flux of ${\boldsymbol{B}} $ due to the Biermann effect is obtainable by solving for the shock part of the "Biermann-only" induction equation, Equation (14). That is to say, we assume the distributional forms of Equations (16), (19), and (26), plug them into Equation (15), and keep only the terms proportional to the Dirac δ. The result may be obtained by inspection of Equation (32), setting the velocities ${\boldsymbol{u}} $ to zero:

Equation (38)

Comparing Equations (38) and (37) and using the continuity of $\nabla {{T}_{e}}\times {\boldsymbol{n}} $, we obtain for the Biermann effect flux of ${\boldsymbol{B}} $ along ${\boldsymbol{n}} $

Equation (39)

This expression is manifestly well defined at discontinuities. As in Equation (32), from which Equation (39) is derived, the only residue of the previously toxic gradient $\nabla {{P}_{e}}$ is contained in the benign direction vector ${\boldsymbol{n}} $, which is effectively $-\nabla {{P}_{e}}/|\nabla {{P}_{e}}|$ at the shock. The effect of the cross product ${\boldsymbol{n}} \times \nabla {{T}_{e}}$ is to eliminate the normal component of $\nabla {{T}_{e}}$, while leaving only the tangential components. It is clear from the form of Equation (39) that only the tangential components of ${\boldsymbol{B}} $ are subject to change according to the R–H condition, as expected.

It is convenient for the purposes of algorithmic implementation to re-express Equation (39) in tensor form, as an antisymmetric tensor ${{{\boldsymbol{G}} }^{(B)}}$ whose components $G_{ik}^{(B)}$ express the flux in the coordinate direction k of the magnetic field component Bi:

Equation (40)

where ${{\epsilon }_{ijk}}$ is the totally antisymmetric Levi–Civita tensor. $G_{ik}^{(B)}$ is, of course, the spatial part of the dual Maxwell tensor due to the Biermann effect. It is worth pointing out that this expression for the flux tensor differs from the "naive" flux

Equation (41)

obtainable directly from the Biermann electric field (Equation (14)) by a pure gradient, which has no effect on the induction equation. It follows that Equation (39) yields a general flux of ${\boldsymbol{B}} $ in the direction ${\boldsymbol{n}} $ that can be used anywhere in the fluid, not only at discontinuities.

The flux in Equation (39) is not the only correction that must be applied to ${\boldsymbol{F}} $. There is also a correction to be applied to the energy part of the flux vector, by virtue of the Poynting flux $(c/4\pi ){{{\boldsymbol{E}} }_{B}}\times {\boldsymbol{B}} $ that arises in connection with the charge-separation electric field ${{{\boldsymbol{E}} }_{B}}$. This term is a bit puzzling at first: since ${{{\boldsymbol{E}} }_{B}}$ is normal to the discontinuity, the Poynting flux is tangential to the discontinuity, and it is not immediately obvious how such a flux is to be pressed into service in an R–H relation.

Nonetheless, the required flux may be inferred in a manner analogous to the calculation by which we derived ${{{\boldsymbol{F}} }_{{\boldsymbol{B}} }}$, above: We solve the "Biermann-only" energy equation

Equation (42)

inserting the distributional forms of Equations (16), (19), (26), as well as

Equation (43)

After collecting terms proportional to the Dirac δ, we obtain

Equation (44)

where in Equation (44) we distinguish between ${\boldsymbol{B}} $—the magnetic field to the left or right of the interface—and $\lt {\boldsymbol{B}} \gt \equiv \frac{1}{2}({{{\boldsymbol{B}} }_{u}}+{{{\boldsymbol{B}} }_{d}})$, the average of the upstream and downstream fields at the interface. Comparing Equation (44) and Equation (37), we obtain for the energy flux along ${\boldsymbol{n}} $ due to the Biermann effect

Equation (45)

Here we see a potential problem: the third term in Equation (45) requires an average of the upstream and downstream currents $\nabla \times {\boldsymbol{B}} $. This is a difficulty if a true flux is to be extracted and pressed into service in a numerical scheme, since possible discontinuities in ${\boldsymbol{B}} $ make such a term ambiguous. The resolution of the difficulty is that, in general, there is always some finite resistivity in real plasmas. As we will see below in Section 4.1, the presence of finite resistivity removes the discontinuity in ${\boldsymbol{B}} ,$ allowing us to replace $\lt {\boldsymbol{B}} \gt $ with ${\boldsymbol{B}} $. We then obtain for the flux of energy in the direction ${\boldsymbol{n}} $

Equation (46)

To summarize, at a cell interface with a normal vector ${\boldsymbol{n}} $, the hydrodynamic fluxes are to be adapted to the Biermann effect by adding to the implemented flux vector ${\boldsymbol{F}} $ an additional flux vector ${{{\boldsymbol{F}} }^{(B)}}({\boldsymbol{n}} )$, given by

Equation (47)

We expect a discretization of the Biermann effect based on these flux expressions to be convergent so long as the mesh resolves the scales on which Te and ${\boldsymbol{B}} $ are continuous. The scale at which Te is continuous is the scale length ${{\lambda }_{T}}$ of the electron thermal conduction precursor region (Shafranov 1957; Jaffrin & Probstein 1964; Zel'dovich & Raizer 2002), discussed in Section 2.2. This length scale is characteristic of the variation of temperature perpendicular to the shock; however, it is also the length scale over which heat diffuses transversely during the time required for the shock to travel a distance ${{\lambda }_{T}}$ (and hence traverse its own thermal precursor). It is therefore the scale to be resolved in order for the discretization of Equation (39) to converge. At coarser resolutions, Te appears as discontinuous at shocks as all the other fluid variables, and the discretization discussed here will not yield converged results for ${\boldsymbol{B}} $. Only as the thermal precursor zone is resolved can convergence be expected.

In order for the discretization to produce convergent results for energy, it is necessary that the simulation resolve the resistive length scale discussed in Section 4.1. Note, however, that for many physical situations of interest, the magnetic field intensities generated by the Biermann effect are so tiny that the correction due to the Biermann effect to the energy flux may justifiably be neglected—it is dwarfed by the fluxes of thermal and kinetic energy, and possibly also by the advected flux of existing magnetic field energy. Under such circumstances, it is an acceptable approximation to neglect the Biermann energy flux, if a non-resistive calculation is of interest, or if the resistive scale cannot be resolved.

4. TWO NOVEL PHYSICAL EFFECTS ARISING FROM SHOCK BIERMANN

We now exhibit two novel predictions of the theory of the Biermann effect at shocks. They are the existence of two magnetic precursors to the shock wave, which lead the wave in regions whose extent depends on the upstream conditions: a resistive magnetic precursor, which arises due to magnetogeneration in the shock "leaking" into the upstream fluid in consequence of the presence of finite nonzero resistivity, and a thermal magnetic precursor, due to smooth-flow Biermann-effect magnetogeneration in the fluid motions set up by the electron thermal precursor. We discuss these in turn.

4.1. The Resistive Magnetic Precursor

When the resistivity η in the induction equation, Equation (3), is finite, a new effect appears in the vicinity of a discontinuity: a resistive magnetic precursor travels ahead of the discontinuity, as the impulsively generated field due to the Biermann effect at the shock diffuses out. This effect is analogous to the well-known electron thermal precursor that precedes a plasma shock, which was discussed in Section 2.2. The structure of the thermal precursor can be estimated by balancing diffusion and advection of thermal energy in the frame of the shock, in the presence of impulsive shock heating (Zel'dovich & Raizer 2002). Similarly, as we will presently see, the structure of the magnetic precursor is a consequence of the balance between magnetic field advection and diffusion in the frame of the shock, in the presence of impulsive magnetogeneration at the shock.

We consider a small portion of the shock surface, which we will treat as approximately planar, and we will work in the rest frame of the shock. We assume an approximately steady state in the shock frame, so that we may set $\partial /\partial t=0$ in the field equations. The induction equation then becomes

Equation (48)

where ${{{\boldsymbol{u}} }^{\prime }}$ is the fluid velocity in the rest frame of the shock.

We will keep only the shock part of the Biermann flux in Equation (48), that is, the right-hand side of Equation (24). In effect, this amounts to assuming that the magnetic field generation from the Biermann effect is impulsive at the shock, so that the smooth-flow Biermann effect is much smaller than the effect at the shock. This approximation is justified if the size ${{\lambda }_{T}}$ of the electron thermal precursor is much smaller than the size of ${{\lambda }_{B}}$ of the magnetic precursor, so that thermal gradients are unimportant except at the discontinuity. The plausibility of this condition will be verified in Section 4.3.

We may adopt the local simplification ${\Psi }({\boldsymbol{x}} ,t)={\boldsymbol{n}} \cdot {\boldsymbol{x}} $, $|\nabla {\Psi }|=|{\boldsymbol{n}} |=1$ for the discontinuity-tracing level function. We further adopt coordinates such that x1 is along ${\boldsymbol{n}} $, so that ${\boldsymbol{n}} ={{{\boldsymbol{e}} }_{1}}$, and such that x2 is along $\nabla {{T}_{e}}\times {\boldsymbol{n}} $. Away from the discontinuity, we ignore all derivatives except for $\partial /\partial {{x}_{1}}$.

The presence of the resistive term changes the discontinuous structure of the solution. This term is the divergence of the resistive flux of ${\boldsymbol{B}} $. That flux behaves analogously to the Fick's law heat flux $-\kappa \nabla T$, in that it opposes gradients in ${\boldsymbol{B}} $. A discontinuous ${\boldsymbol{B}} $ is disallowed, because it would result in an infinite restoring resistive flux. ${\boldsymbol{B}} $ is therefore now continuous. By the R–H condition for transverse momentum (Gurnett & Bhattacharjee 2005, Chapter 7), it follows that ${\boldsymbol{u}} _{T}^{\prime }$ is also continuous. Only $u_{n}^{\prime }$ and Pe have solution discontinuities in this case. Assuming that the upstream fluid is at rest in the lab frame, we therefore set ${{{\boldsymbol{u}} }^{\prime }}=u_{n}^{\prime }{{{\boldsymbol{e}} }_{1}}$, with

Equation (49)

By the coordinate choice, assuming that the far upstream fluid is unmagnetized, and by virtue of Equation (48), we may set ${\boldsymbol{B}} =B{{{\boldsymbol{e}} }_{2}}$. Putting this all together, we obtain the following structure equation for the magnetic precursor:

Equation (50)

The meaning of this equation is that the precursor structure is determined by the balance of magnetic diffusion and magnetic advection, in the presence of impulsive magnetic field generation due to the Biermann effect.

In the upstream (${{x}_{1}}\gt 0$) region, we may neglect gradients in η (which depends on Te). Equation (50) becomes

Equation (51)

the solution of which, given ${{{\rm lim} }_{{{x}_{1}}\to +\infty }}B({{x}_{1}})=0,$ is

Equation (52)

Here B0 is an integration constant. We see that the precursor has an exponential shape and a characteristic length ${{\lambda }_{B}}$ given by

Equation (53)

Here $\eta =\eta ({{T}_{e,+\infty }})$ is the value of the resistivity far upstream of the discontinuity.

We may obtain a relation for B0 from the jump condition at the discontinuity, by integrating Equation (50) across a vanishingly small shock crossing path. The result is

Equation (54)

where now $\eta ({{T}_{e}})$ is evaluated at the shock. This condition, together with the upstream solution and boundary condition, determines B0.

The ideal MHD jump conditions may be recovered in the limit of spatially constant η with $\eta \to 0$. In that case, the downstream solution satisfying finite-B boundary conditions at ${{x}_{1}}\to -\infty $ is ${{B}_{d}}({{x}_{1}})={{B}_{0}}$, and the first term in Equation (54) is just $-D{{B}_{0}}$. Equation (33) follows immediately.

4.2. The Electron Thermal Precursor

One further consequence of the presence of the thermal precursor in Te described in Section 2.2 is that, in general, there will be some small amount of magnetic field generated by the smooth-flow Biermann effect ahead of the shock, owing to plasma motions in the preheated region, whose size is ${{\lambda }_{T}}$, given by Equation (6). In general, the field intensity due to this effect can be expected to be small compared to the field intensity due to the resistive magnetic precursor, since the gradients in the thermal magnetic precursor are small compared to those available near the shock itself. In the constant-conductivity simulations that we discuss in Section 6.2.4, we will see that the field intensity is in fact quite small.

Nonetheless, this smallness does not necessarily preclude the possibility that the thermal magnetic precursor might be experimentally observable. Thermal conductivity upstream of a plasma shock is, in general, not constant, but rather depends strongly on electron temperature—${{\kappa }_{e}}\sim T_{e}^{5/2}$ (Spitzer 1962). The actual structure of the thermal precursor is not the gentle exponential decay that one expects for a constant ${{\kappa }_{e}}$, but rather exhibits the steep gradient near its outer terminus shown in Figure 7.20 of Zel'dovich & Raizer (2002). It is possible that this large gradient may be responsible for an enhancement of the Biermann effect at the terminus of the thermal conduction precursor zone. This interesting possibility is beyond the scope of this paper.

One consequence of the presence of a thermal magnetic precursor is that even in a nonresistive approximation, the notion of an unmagnetized fluid upstream of the shock, on which, for example, Equation (33) is based, is not strictly correct. It is a valid approximation only when the Biermann generation rate in the thermal precursor is small compared to that of the shock.

4.3. Physical Length Scales

It is useful to establish some of the relevant length scales under conditions of interest. Below, we establish these scales for conditions prevailing in HEDP experiments. In Section 5 we will establish the scales characteristic of galaxy cluster formation.

Using the expression for η from Huba (2007) in the expression for the resistive magnetic precursor length ${{\lambda }_{B}}$ given in Equation (53), we have the following expression for ${{\lambda }_{B}}$ in a fully ionized plasma:

Equation (55)

where Λ is the usual Coulomb logarithm. Similarly, the expression for the electron thermal precursor length scale ${{\lambda }_{T}}$ is (Zel'dovich & Raizer 2002, Chapter 7, Section 12)

Equation (56)

where ${{T}_{e,0}}$ is the electron temperature at the shock, and where ξ is a number in the range 1–2.

The numbers in Equations (55) and (56) have been scaled to conditions that are routinely obtainable in large laser facilities such as Omega and NIF (see, e.g., Gregori et al. 2012). In these conditions, it is clear that the assumption ${{\lambda }_{T}}\ll {{\lambda }_{B}}$, required for the validity of the derivation of ${{\lambda }_{B}}$, is easily satisfied.

It is also immediately apparent that the magnetic precursor length scale is a macroscopic scale that can, in principle, be well tailored to the physical dimensions of the interaction chambers of such facilities. A carefully designed experiment, which launches a shock into a relatively cold upstream plasma (note the dependence of ${{\lambda }_{B}}$ on Te) should in principle be able to detect the resistive precursor due to the Biermann effect.

Note that as the value of Te upstream rises, ${{\lambda }_{B}}$ decreases as $T_{e}^{-3/2}$, while ${{\lambda }_{T}}$ increases as $T_{e}^{5/2}$. We therefore only need to increase the upstream temperature to about 7 eV for ${{\lambda }_{T}}$ to become comparable to ${{\lambda }_{B}}$, and at warmer temperatures than this ${{\lambda }_{T}}$ becomes dominant. It is therefore plausible that a physical situation could be created in which the thermal magnetic precursor might be observable without interference from the resistive magnetic precursor.

5. THE BIERMANN EFFECT AND GALAXY CLUSTER FORMATION

In this section we attempt to estimate the strength of the seed magnetic fields that result from a correct treatment of the Biermann effect at shocks in galaxy clusters at the time of magnetogenesis ($z\sim 5-3$) and the physical length scales of the resistive magnetic precursor and the electron thermal precursor in this context.

5.1. Protogalaxy Cluster Field Generation

As discussed in the Introduction, the Biermann effect has been invoked as the source of seed magnetic fields that can be amplified by the turbulent dynamo mechanism (Kulsrud et al. 1997).

Given the fact that in these studies the gradients near shocks that were used to calculate the Biermann effect field generation rates are artifacts of the hydrodynamic advance schemes, it is possible—even likely—that the calculated field strengths are incorrect. We now attempt to estimate the typical strength of magnetic fields that are to be expected in early galaxy cluster formation, based on a corrected treatment of the Biermann effect. We do not perform new simulations, but rather attempt to estimate the typical field strength based on published simulation data, in a preliminary effort to determine the reliability of field strengths inferred to date in studies of galaxy cluster formation.

While this is not an easy task without actual simulation data on hand to analyze, the analysis work described in Miniati et al. (2000), based partly on ΛCDM simulations described in Cen & Ostriker (1994), provides enough information for us to make a rough estimate of the typical field strength. The simulation in question has the properties ${{H}_{0}}=60$ km s−1 Mpc−1, ${{{\Omega }}_{M}}=0.45$, ${{{\Omega }}_{{\Lambda }}}=0.55$, ${{{\Omega }}_{b}}=0.043$. We now give a relatively detailed description of our analysis of the information in Miniati et al. (2000), in order to make clear both the basis for the estimated field strength and the considerable uncertainty that attend these estimates.

Our starting point is Equation (33), which gives the postshock field strength due to the Biermann effect, assuming an unmagnetized upstream fluid. In order to use this equation to estimate field strengths, we need typical values of the shock speed D, the compression ratio ${{\rho }_{d}}/{{\rho }_{u}}$, and the tangential gradient $|{{\nabla }_{\bot }}({{k}_{B}}{{T}_{e}})|\equiv |\nabla ({{k}_{B}}{{T}_{e}})\times {\boldsymbol{n}} |$ that arise at primordial epochs.

For the sake of exploiting the information available in Miniati et al. (2000), we will choose z = 3 as our fiducial primordial epoch. Examining the top right panel of Figure 8 in Miniati et al. (2000), we conclude that a not-atypical value of $|{{\nabla }_{\bot }}({{k}_{B}}{{T}_{e}})|$ is

Equation (57)

Neglecting momentarily the fact that the simulation certainly fails to resolve the thermal precursor length scale ${{\lambda }_{T}}$ (Equation (60) below), the uncertainty in this estimate seems to be a factor of 2 or so.

Miniati et al. (2000) supply shock speeds at redshift z = 0 (see their Figure 5). These may be approximately shifted to z = 3 using Figure 10 of Miniati et al. (2000), which shows that the kinetic energy processed by shocks increased by a factor of 15 between z = 3 and z = 0. This suggests that temperatures increased by about that much (neglecting shock filling factor differences between the two redshifts) and that velocities increased by about a factor of 4. We had typical shock temperatures of $\sim {{10}^{6}}$ K at z = 3 in the data leading to the estimate in Equation (57). This corresponds to a temperature $\sim {{10}^{7}}$ K at z = 0. From the top right panel of Figure 5 of Miniati et al. (2000), this corresponds to a shock speed of about $6\times {{10}^{7}}$ cm s−1 at z = 0, and hence $D\sim {{10}^{7}}$ cm s−1 at z = 3.

Putting these numbers in Equation (33) and assuming the strong-shock limit ${{\rho }_{d}}/{{\rho }_{u}}=4$ appropriate to a $\gamma =5/3$ ideal gas, we obtain

Equation (58)

This value should be regarded as a lower limit since, as remarked above, the simulation does not resolve ${{\lambda }_{T}}$. An upper bound on B can be obtained by replacing, in Equation (57), the 5 Mpc length scale estimated from Figure 8 of Miniati et al. (2000) with ${{\lambda }_{T}}$ from Equation (60). This gives $B\lesssim {{10}^{-19}}$ G, i.e., a value that is a factor of about 103 higher. Obviously, this is a very conservative bound, as it assumes temperature fluctuations of order unity over the diffusive scale ${{\lambda }_{T}}$.

It is clear from the very uncertain nature of the factors used to construct this estimate that the field strength given in Equation (58) is itself subject to considerable uncertainty. A few actual MHD simulations with a correct implementation of the Biermann effect seem required to establish how much field strength is in fact made available by the Biermann effect, for turbulent dynamo effects to amplify. Such simulations would be challenging given the spatial scale ${{\lambda }_{T}}$ that needs to be resolved.

5.1.1. Galaxy Cluster Formation Length Scales

We now estimate the physical length scales of the resistive magnetic precursor and the electron thermal precursor for the case of galaxy cluster formation using the physical parameters inferred in Section 5.1 from Miniati et al. (2000) for galaxy cluster formation at $z\simeq 3$. As in Section 5.1, we use ${{T}_{e}}={{10}^{6}}$ K as our reference temperature and $D={{10}^{7}}$ cm s−1 as our reference shock speed. The assumed cosmological parameter ${{{\Omega }}_{b}}=0.043$ corresponds to a proton number density ${{n}_{i}}\approx {{10}^{-7}}$ cm−3 at z = 0, and hence to ${{n}_{i}}\approx 6\times {{10}^{-6}}$ at z = 3, which we take as the density upstream of the shocks. With these parameters, the Coulomb logarithm is ${\rm ln} {\Lambda }\approx 40$. Setting Z = 1, A = 1, we then have

Equation (59)

and

Equation (60)

The much more tenuous and considerably hotter primordial plasma reverses the relative sizes of ${{\lambda }_{B}}$ and ${{\lambda }_{T}}$ compared to the HEDP case considered above: ${{\lambda }_{B}}$ is now negligible, while ${{\lambda }_{T}}$ is an astrophysically significant length. It is worth comparing ${{\lambda }_{T}}$ to ${{\lambda }_{e}}$ and ${{\lambda }_{i}}$, the expected Spitzer mean free paths of electrons and ions in the plasma, as a consistency check. This is the product of the mean collision time τ with the thermal speed ${{(3{{k}_{B}}T/m)}^{1/2}}$. Using the expression for ${{\tau }_{e}}$ from Huba (2007), we have

Equation (61)

The mean free path is independent of particle mass, and for Z = 1 it is the case that ${{\lambda }_{e}}={{\lambda }_{i}}$. This length is the characteristic size of the viscous shock sheath and is reassuringly shorter than the semihydrodynamically scaled ${{\lambda }_{T}}$.

It is also worth considering whether the plasma is collisional, as we have implicitly assumed by using Spitzer-type transport coefficients. The magnetic field strength estimated in Section 5.1 gives rise to electron gyroradii ${{\lambda }_{c}}$ of characteristic size

Equation (62)

Since ${{\lambda }_{c}}\sim 100{{\lambda }_{e}}$, the flow is comfortably collisional for the chosen physical parameters. An increase in the estimate of B by two orders of magnitude or more would change this conclusion; however, it should be noted that at least at the outset of the process of cosmic magnetogenesis envisioned here, field strengths were probably small enough that the collisional plasma approach is valid. Primordial field strengths that may have originated in early-universe phase transitions or in inflationary scenarios are poorly constrained by theory and observation, but they are believed to have plausible values $B\lt {{10}^{-22}}$ G in all but a few scenarios (see Widrow 2002; Widrow et al. 2012). It follows that, initially at least, the Spitzer-type transport coefficients that we have employed here are likely to be appropriate. This is in contrast to the apparent suppression of electron conductivity at later times (see, e.g., Markevitch et al. 2003; Russell et al. 2012; Sanders et al. 2013; ZuHone et al. 2013).

6. NUMERICAL VERIFICATION

6.1. Implementation Using the FLASH Code

We have implemented the corrected Biermann effect algorithm, described in Section 3.4, within the publically available FLASH hydrodynamic simulation framework (Fryxell et al. 2000; Dubey et al. 2009; Tzeferacos et al. 2014). FLASH is a modular and extensible multiphysics scientific simulation software package that has been widely used for reactive compressible flows typical of astrophysical situations, HEDP applications, cosmology, computational fluid dynamics, and fluid–structure interactions. In particular, FLASH makes available both a two-temperature single-fluid model and resistive MHD, which makes it an ideal platform for testing the proposed algorithm.

Below we furnish an outline of the workings of the code, including the newly implemented Biermann effect algorithm. Much more detail about FLASH is available in the FLASH Users Guide.1 We plan to include an implementation of the new Biermann effect algorithm in a future release of the code.

The two-temperature MHD model that we use for our numerical tests is expressed in conservation form by the following dynamical system:

Equation (63)

Equation (64)

Equation (65)

Equation (66)

Equation (67)

These equations differ from Equations (8)–(12) only in that the Biermann term in Equation (65) replaces the one in Equation (10) to reflect Equation (46), while the Biermann term in Equation (67) replaces the one in Equation (12) to reflect Equation (40).

FLASH advances the MHD equations using an unsplit staggered mesh (USM) scheme that prevents the development of magnetic monopoles due to numerical noise in the induction equation, and which is described in Lee (2013).

The thermal conduction and heat-exchange source terms on the right-hand side of Equation (66) are computed in a time-split manner, separately from the ideal MHD advance. In particular, thermal conduction advance is fully implicit and works as described in Section 17.1.4 of the FLASH Users Guide, while heat exchange operates as described in Section 16.5 of the FLASH Users Guide. The remaining, advective part of Equation (66) is implemented by treating se as a passively advected mass scalar. This equation requires an EOS implementation capable of using electron entropy as an input variable. For the current case of perfect gas EOS and total ionization, the standard Sackur–Tetrode equation for entropy is implemented in the EOS.

The resistive terms in Equations (65) and (67) are not treated implicitly, but rather are added directly to the MHD fluxes and are thus treated explicitly, placing diffusive stability limitations on the maximum time step. The Biermann terms in Equations (65) and (67) are also added explicitly to the MHD fluxes. These flux modifications are performed in such a way as to preserve the solenoidal character of the magnetic field, by adding them to the ideal MHD Godunov fluxes before these are interpolated to the edge-centered electric fields of the USM scheme, after which each face-centered normal magnetic field component is updated by the circulation integral of the electric field along the edges bounding the face (Lee 2013).

We have not performed an analysis of the time step limitation imposed by the Biermann effect, trusting instead that the small magnitude of the effect in the cases considered makes it unlikely that Biermann time step constraints could dominate those due to the hyperbolic and diffusive terms. In this connection, it is noteworthy that the Biermann term, being quadratic in fluid derivative terms, makes no contribution to wave dispersion relations obtained by linearization about uniform solutions and therefore does not appear to affect ordinary wave speeds. Intuitively, one would therefore expect any time step limitations due to the Biermann effect to be higher-order corrections, hardly affecting the numerical evolution of plasmas for which the Biermann term is small in magnitude. Nonetheless, a full analysis of the effect of the Biermann term on the hyperbolic structure of the system of PDEs would be useful—possibly even necessary—for cases of intense field generation. Such an analysis would be somewhat complicated by the non-quasi-linear structure of the Biermann term, necessitating an extension of the PDE system to quasi-linear normal form.

6.2. The Simulations

The simulations presented below represent idealized situations simplified for the sake of verifying the code and for illustrating the numerical and physical principles under discussion. In particular, the conductivity ${{\kappa }_{e}}$, electron–ion equilibration timescale ${{\tau }_{ei}}$, and resistivity η are all treated as constants. FLASH does have the capability to use Spitzer-type functions of the thermal state for these parameters, but this would complicate the results unnecessarily without adding any real value to these verification tests.

In what follows, we assume fully ionized hydrogen—A = 1, Z = 1, adiabatic index $\gamma =5/3$, ${{c}_{v,e}}=\frac{3}{2}{{k}_{B}}{{N}_{A}}$, where NA is Avogadro's number and kB is Boltzmann's constant. Simulations are conducted in cylindrical coordinates, assuming azimuthal symmetry, and are therefore two-dimensional. In addition, we impose a reflection boundary on the R-axis (where R is the perpendicular distance from the axis of rotational symmetry, which is to say, R is the cylindrical radius) so that the domain represents a hemisphere of a solution with reflection symmetry about that axis. In every case, the domain is a 4 cm radius cylinder that extends 4 cm in the z direction from the R axis. The boundary conditions at R = 4 cm and at z = 4 cm are outflow.

All of the verification tests that we present are variations on the theme of a Sedov-esque explosion: an analytic self-similar Sedov solution is (1) modified to prevent temperatures from diverging at the center by setting all variables constant inside some chosen radius; (2) smoothed with a Gaussian near the shock, to ease instabilities that can otherwise appear near the shock; (3) where required, distorted to an ellipsoidal profile capable of producing real Biermann effect fields; and (4) scaled to velocities and pressures estimated to produce usable simulation times given the resolutions and domain size in play. Values of shock ellipticity, ${{\kappa }_{e}},$and η are then chosen to produce reasonable thermal and magnetic precursor region sizes and reasonable field generation rates, and values of ${{\tau }_{ei}}$ are chosen to be harmless. In every case, we run an initial model for a time without the Biermann effect, so as to allow transients in the two-temperature hydrodynamic variables to subside. Then we restart the calculation with Biermann field generation turned on.

Magnetic fields are always azimuthal in these verification tests—since the gradients of Te and Pe are always in the $R-z$ plane, the Biermann effect only generates a nonzero field along the azimuthal direction.

6.2.1. Null Test: A Spherical Shock

It is not a simple matter to construct a nontrivial analytic verification solution of a plasma generating magnetic field through the Biermann effect. A trivial solution, on the other hand, may be straightforwardly constructed by imposing symmetry requirements that align the gradients $\nabla {{P}_{e}}$ and $\nabla {{T}_{e}}$, thus guaranteeing zero-field generation. We use a spherical Sedov-like explosion as an example of such a test. This is in fact the test that revealed the Biermann catastrophe in the first place (Fatenejad et al. 2013).

We assume an initial shock radius of 1.15 cm, in an ambient medium with $\rho =1$ g cm−3 and P set to a negligible value. We choose an initial velocity scale inside the shock that leads to a shock velocity that is on average about $3.2\times {{10}^{4}}$ cm s−1 during the course of the simulation. Other physical parameters are ${{\kappa }_{e}}=1.0\times {{10}^{12}}$ erg K −1 cm−1 s−1 and ${{\tau }_{ei}}=2\times {{10}^{-14}}$ s. The initial solution is advanced for 10−5 s with no field generation and then restarted and run for another 10−5 s with field generation turned on. We repeated this experiment with the correct Biermann flux term of Equation (40), with the naive flux of Equation (41), and with the time-split source term of Equation (15), so as to compare the performance of the three algorithms.

The results of this study are displayed in Figure 2. The top left panel is a color map view of the electron temperature Te. The dashed shock-crossing line across the bottom of this figure is the location of the data displayed in the top right panel, which shows electron and ion temperatures. The shock is visible as the peak in Ti. The thermal precursor is clearly visible, as is the continuous behavior of Te, which changes slope at the shock in accordance with standard plasma shock theory (Zel'dovich & Raizer 2002, Chapter 7, Section 12). The ion temperature also rises in the precursor region in consequence of the electron–ion equilibration term.

Figure 2.

Figure 2. Two-dimensional simulation of a spherical shock in a two-temperature hydrogen plasma. The images refer to the final time step, at $t=2\times {{10}^{-5}}$ s, and to the highest-resolution simulation. Top left: electron temperature distribution. Top right: electron and ion temperatures along the shock-crossing line shown at the bottom of the figure in the top left panel. Bottom left: ${{B}_{\phi }}$ due to numerical noise. Bottom right: total magnetic energy as a function of simulation resolution.

Standard image High-resolution image

The bottom left panel of Figure 2 shows the magnetic field strength generated by the correct treatment of Equation (40). There is a spurious nonzero field due to numerical noise that is clearly being generated. However, the bottom right panel shows that this field is converging to zero with increasing resolution. This figure displays the square root of the total magnetic energy in the domain. Since the analytic solution for the magnetic field strength is zero everywhere, this quantity functions as an un-normalized L2-norm of the difference between the analytic and numerical solutions. The blue circles show the convergence with resolution of the correct formulation. This convergence behavior is in contrast to the behavior of the other two algorithms, which fail altogether to converge.

6.2.2. Verification with an Ellipsoidal Shock

Next we exhibit simulations designed to produce nonzero values of the magnetic field strength. We start with an ellipsoidal shock surface, obtained by distorting the spherical Sedov solution. The initial configuration has a semimajor axis (aligned with the R-axis) of 2 cm and a semiminor axis (aligned with the z-axis) of 1.667 cm. The ambient medium is uniform with $\rho =1$ g cm−3 and P set to a negligible value. We choose a velocity scale that leads to shock velocities of about $6\times {{10}^{6}}$ cm s−1. The higher shock speed is chosen to produce somewhat intense magnetic field strengths. To keep the thermal precursor resolved, we increase the conductivity to ${{\kappa }_{e}}=1.0\times {{10}^{13}}$ erg K−1 cm −1 s−1. We also set ${{\tau }_{ei}}={{10}^{-14}}$ s. The initial solution is advanced for $4.86\times {{10}^{-8}}$ s with no field generation and then restarted and advanced to a simulation time $t=3\times {{10}^{-7}}$ s with field generation turned on.

The advance of the shock during the period when magnetic field is generated is illustrated by the two pressure color map plots in the top panels of Figure 3. The bottom left panel of Figure 3 shows the magnetic field intensity generated according to the correct flux of Equation (40), ranging into the tens of gauss. The solid red line in the figure shows the shock location, while the red arrows display the shock unit normal vector. The bottom right panel displays the difference between the correct and the naive flux formulations, which is substantial at the outgoing shock surface.

Figure 3.

Figure 3. Ellipsoidal shock. Top left: initial pressure distribution. Top right: final pressure distribution. Bottom left: final magnetic field distribution. The solid red line displays the location of the shock, while the red arrows display the direction of the shock normal. Bottom right: difference between final magnetic field distributions computed using the correct and incorrect Biermann flux.

Standard image High-resolution image

We do not have an analytical solution for the magnetic field distribution to compare to the output of these simulations. We do, however, have the relation between shock quantities expressed by Equation (32), which determines the jump condition for ${\boldsymbol{B}} $. By locating points adjoining the shock and computing the local shock velocity at those points, we can then verify that Equation (32) in fact obtains, to some accuracy limited by the numerical approximation. This is, in principle, a nontrivial verification, since the code does not know about the jump condition Equation (32) as such—it only knows the fluxes expressed by Equation (40), which were inferred from the jump condition. Recovery of the jump condition thus constitutes a verification that the code correctly models the theory.

To perform this verification, we first locate cells adjoining the shock by the method described in the Appendix, which yields a one-cell-wide shock surface by fitting the mass, momentum, and energy R–H conditions to the neighboring data, using a speed-of-sound weighted inner product on the state space, and treating the shock speed D as a free fit parameter. The shocked cells are the ones that fit the R–H conditions with small fit residuals (${{R}^{2}}\lt 0.01$), together with some other natural auxiliary consistency conditions described in the Appendix. The fitted shock speed is then used in the verification of Equation (32).

In the absence of normal component field Bn, Equation (32) becomes

Equation (68)

where we have defined "advection" terms ${{a}_{u,d}}$ and "Biermann" terms ${{b}_{u,d}}$ by

Equation (69)

Equation (70)

Equation (71)

Equation (72)

The sum of terms in Equation (68) is required to be zero. In a discretized numerical PDE integration, this really means that the terms must cancel up to some truncation or rounding precision, which is expressed relative to the largest of the magnitudes in play. We therefore define the "shock condition parameter" C by

Equation (73)

We calculate the value of C at cells along the shock front, at each of our four resolution levels, for both the correct flux and the "naive" flux implementations of the Biermann effect. We display cumulative distributions of C in the top panels of Figure 4. It is evident from these figures that the distribution of C is in fact centered near zero for both flux implementations. The widths of the distributions behave very differently, however. In the case of the correct Biermann flux implementation, the distributions get narrower with each refinement of resolution, providing some evidence of convergence with resolution to the expected result. In the case of the naive flux implementation, there is no such evidence of convergence. The same observations can be made of the plots in the middle panels of Figure 4, which summarize the means and sample standard deviations of the C-distributions, as a function of resolution. The convergence of the correct flux implementation and the convergence failure of the naive implementation are clear in these figures.

Figure 4.

Figure 4. Left: correct Biermann flux term. Right: "naive" Biermann flux term. Top: cumulative distribution of the normalized magnetic shock condition C, defined by Equation (73), evaluated at points along the shock surface, for four different resolutions, illustrating convergence to the correct jump condition. Middle: sample standard deviations of C, as a function of resolution. Bottom: total magnetic energy as a function of simulation time for the four resolutions studied.

Standard image High-resolution image

The bottom panels of Figure 4 show the evolution of total magnetic energy in the domain as a function of time, for the different resolutions and for the two flux implementations. The correct flux implementation appears to be converging (although not in any strong sense converged) by this measure, even at late time, whereas any evidence of convergence in magnetic field energy is simply lacking for the naive flux implementation.

6.2.3. The Resistive Magnetic Precursor

In this set of simulations, we maintain the simulation parameters described in Section 6.2.2 but also turn on the resistivity η to a finite positive value, $\eta =7.8\times {{10}^{5}}\;{{{\rm s}}^{-1}}$, chosen in conjunction with the other parameters to yield an easily discernible exponentially decaying resistive precursor described by Equation (52). We repeat the simulation strategy of Section 6.2.2, letting the simulation advance to $t=4.86\times {{10}^{-8}}$ s with no field generation, and then restarting it and advancing to a simulation time $t=3\times {{10}^{-7}}$ s with field generation (using the correct flux formulation) turned on. This advance time is adequate for the establishment of the expected precursor structure, as we demonstrate below.

The left panel of Figure 5 shows a color map of the distribution of magnetic field strength across the domain. The shock location is shown by the solid red line, and the superposed vectors indicate the shock-normal direction. The magnetic field is evidently smoothed out by the resistivity, as can be seen by a direct comparison with the bottom left panel of Figure 3. The precursor is already evident in this figure, as the substantial amount of magnetic field that has "leaked" ahead of the shock.

The diagonal white dashed shock-crossing line in the left panel of Figure 5 illustrates the line along which data were extracted to produce the right panel of Figure 5. This figure shows the magnetic field values plotted along that line. The vertical solid red line at X = 0 marks the location of the shock. The exponential decay of the field is easily visible in this figure.

At each shock location, the value of η can be combined with the local shock velocity to calculate ${{\lambda }_{{\rm Predicted}}}$, the expected decay length of the precursor, according to Equation (53). At the same time, the actual decay length ${{\lambda }_{{\rm Measured}}}$ can be measured directly by comparing the field strength at two suitably selected distances along the local shock normal. We have performed both calculations at each shock location and compared them. The results are displayed in Figure 6, where we plot the cumulative distribution of ${{\lambda }_{{\rm Measured}}}/{{\lambda }_{{\rm Predicted}}}$. The distribution appears to be consistent with a mean value of 1, as expected, with some scatter. The scatter is not surprising, since the expression in Equation (53) for ${{\lambda }_{{\rm Predicted}}}$ is an approximation, assuming, as it does, propagation of magnetic field from a planar shock. Since the shock is, of course, not planar, the value of the Biermann generation rate changes appreciably across the shock over distances not hugely different from ${{\lambda }_{{\rm Predicted}}}$ itself. Under the circumstances, then, the level of agreement between observation and prediction is satisfactory.

Figure 5.

Figure 5. Left: magnetic field distribution due to passage of shock, in the presence of finite resistivity. The simulation differs from the one shown in Figure 3 only by the presence of nonzero resistivity. The solid red line shows the location of the shock, while the vectors represent the unit normal to the shock. The magnetic precursor can be clearly seen ahead of the shock surface. The white dashed shock-crossing line illustrates the location of the data displayed in the right panel. Right: magnetic field strength magnitude along the white dashed shock-crossing line shown in the left panel. The position of the shock is shown by the vertical solid red line. The predicted exponential decay of the precursor is evident.

Standard image High-resolution image
Figure 6.

Figure 6. Cumulative distribution of the measured magnetic precursor decay length, normalized to the predicted length, at 849 points along the shock surface. The mean is shown by the dotted line, while the standard deviation is illustrated by the error bar.

Standard image High-resolution image

6.2.4. The Thermal Magnetic Precursor

In our final simulation, we explore the properties of the thermal magnetic precursor by returning to a nonresistive simulation, similar to the ones discussed in Section 6.2.2, and differing from them only in that ${{\kappa }_{e}}$ is 10 times larger: ${{\kappa }_{e}}=1.0\times {{10}^{14}}$ erg K−1 cm−1 s−1. This broadens the thermal precursor zone to about 0.2 cm, making it easier to discern.

The results at the final time step are illustrated in Figure 7. The top left panel shows the distribution of magnetic field in the domain. The solid red line in the figure shows the location of the shock, while the arrows show the direction of the shock normal. In the top right panel we see the electron and ion temperatures along the shock-crossing line shown in the previous panel. The shock location coincides with the sharp drop in Ti, while the thermal conduction precursor zone, which has an exponential structure for the constant conductivity model used here, may be clearly seen in the upstream behavior of Te.

Figure 7.

Figure 7. Ellipsoidal shock, with 10 times greater thermal conductivity. Top left: final magnetic field distribution. The solid red line shows the shock location, while the arrows show the direction of the shock normal. The white dashed shock-crossing line illustrates the location of the data displayed in the remaining panels. Top right: electron and ion temperatures along the white dashed shock-crossing line in the previous panel. The shock location and thermal precursor are evident. Bottom left: B along the shock-crossing line. Bottom right: semilog plot of $|B|$ along the shock-crossing line.

Standard image High-resolution image

The magnetic field structure along the shock-crossing line is shown in the bottom left panel of the figure. The precursor is difficult to see in this plot, so the bottom right panel shows $|B|$ on a semi-log plot. The magnetic precursor structure is seen in this figure. It clearly has a much lower field intensity than the fields generated by the Biermann effect at the shock. We emphasize again, however, that in part this is due to the highly simplified constant-conductivity model adopted for this work. A true Spitzer-type conductivity, with a $T_{e}^{5/2}$ dependence, would create much sharper gradients at the upstream termination of the precursor zone, potentially generating more intense fields—relative to those generated at the shock—than those shown here.

7. SUMMARY AND DISCUSSION

To summarize our findings: using the kinetic theory of plasma shocks, we have given a description of the Biermann effect within ion viscous shocks. Using this description, we have shown that the convergence failures in the presence of shocks of MHD codes implementing the Biermann effect are not traceable either to a mathematical mis-specification of the Biermann term—which is well defined despite appearances—or to a failure of the Biermann MHD source term to correctly approximate the expected physics. The failure is instead due to naive discretization, a condition that we have shown is curable by exploiting the continuity of the electron temperature Te.

We have described a convergent algorithm that incorporates the Biermann effect within shocks in numerical MHD and implemented it in the FLASH code. We have developed and used verification tests to demonstrate formal convergence using a null solution from spherical shock and to verify predictions for physical conditions near a shock using an ellipsoidal shock. Comparisons of the new algorithm, which provides a correct and accurate treatment of the Biermann effect within shocks, with the previous, naive one exhibit clearly desirable physical and numerical properties present in the new algorithm that were previously wanting.

We have noticed, described, and simulated two previously unrecognized physical effects: a resistive magnetic precursor that leads the shock in the presence of nonzero resistivity, and that is analogous to the well-known thermal precursor caused by the presence of finite electron thermal conduction; and a thermal magnetic precursor, produced by plasma motions in the electron thermal conduction precursor. We have estimated the expected magnitude of both effects in conditions encountered in laboratory experiments at laser facilities and shown that the characteristic length of the two effects can be made macroscopic. In particular, the resistive magnetic precursor is physically measurable in an experimental setup that measures both the position of the shock front and a time series of the intensity of the magnetic field at some location initially ahead of the shock. Such an experiment would have to somehow ensure that the upstream fluid remains unheated by shock-generated radiation and by laser absorption, so as to keep the value of the resistivity high enough to sustain a macroscopically scaled precursor length. The thermal magnetic precursor may be more challenging to observe, as it is intrinsically weaker than the resistive precursor (relative to the field intensity generated at the shock). Further studies with Spitzer-type electron thermal conductivity are needed to determine whether sharp gradients at the upstream termination of the electron thermal precursor can give rise to sufficiently intense magnetic field strengths to be experimentally measurable.

We emphasize again that the requirement for convergence in ${\boldsymbol{B}} $ of the new Biermann effect algorithm is that the electron thermal precursor of the shock should be resolved, while the requirement for convergence in magnetic energy of the new algorithm is that the resistive precursor should be resolved. Of these two requirements, the first is more essential, because neglect of the contribution of the Biermann effect to the flux of magnetic energy is frequently a valid approximation. Demanding resolution of the electron thermal precursor is not a trivial requirement in realistic simulations. It certainly mandates some kind of Adaptive Mesh Refinement (AMR) strategy to resolve the region near the shock. But irrespective of AMR-related efficiencies, merely resolving that length scale can constrain the code to taking very short time steps owing to the Courant–Friedrichs–Lewy stability condition, unless an implicit advancement scheme is implemented for MHD. This issue is the subject of current study, but it lies beyond the scope of the present paper.

With the discussion of the modification due to the Biermann effect of the R–H jump condition on B—Equation (32) and text thereabout—it now seems worth returning to the question of vorticity ${\boldsymbol{ \omega }} $, and the beguiling Equation (5) that suggests its proportionality to ${\boldsymbol{B}} $. Recall that much hangs on the jump conditions—if the jump in ${\boldsymbol{ \omega }} $ maintains the same proportionality to the jump in ${\boldsymbol{B}} $ as do their respective evolution equations, then, subject to some cavils about boundary conditions, ${\boldsymbol{B}} $ is simply ${\boldsymbol{ \omega }} $, at least until resistivity and/or viscosity manifest themselves, and in effect, unmagnetized gasdynamics contains within it the magnetic degrees of freedom of MHD, encoded in derivatives of the velocity field.

It can be seen more clearly now why this implausible assertion is in fact false. In the first place, ${\boldsymbol{B}} $ and ${\boldsymbol{ \omega }} $ are rather different types of quantities from the point of view of kinetic theory, in that ${\boldsymbol{B}} $ is perfectly well defined at kinetic scales—such as inside an ion viscous shock sheath—whereas ${\boldsymbol{ \omega }} $ cannot even be defined consistently in regions where the fluid picture breaks down. The jump condition on ${\boldsymbol{B}} $ is a straightforward consequence of the induction equation and is essentially kinematic in origin, as are the other R–H conditions. The jump condition on vorticity, on the other hand, cannot be inferred from the dynamical PDEs alone, but depends in an essential manner on the EOS (Kevlahan 1997). It would be possible, in principle, to infer a "flux" of ${\boldsymbol{ \omega }} $ from Equation (4) and to construct a "jump condition" using that flux in the R–H conditions, Equation (37). That jump condition would disagree with the one obtained in Kevlahan (1997), as it is manifestly independent of the EOS. The EOS dependence of the jump condition on vorticity implies that even if there exists an EOS according to which the jump condition on ${\boldsymbol{ \omega }} $ preserves the required proportionality to the Biermann-laden jump in ${\boldsymbol{B}} ,$ that proportionality would certainly not be preserved for any other EOS.

In general, then, the passage of a shock certainly spoils Equation (5). This should be no surprise, as the manipulations of the hydrodynamic equations of motion required to arrive at Equation (4) constitute precisely the sorts of prestidigitation with vector derivatives that break down at fluid discontinuities (see, e.g., LeVeque 2002). While Equation (4) may be derived from the momentum equation on either side of the shock, the connection between ${\boldsymbol{B}} $ and ${\boldsymbol{ \omega }} $ is certain to be lost at the shock itself, well before resistive and viscous effects can assert themselves. We conclude that the vorticity connection is not a useful tool for studying magnetogeneration by the Biermann effect in the presence of shocks.

In concluding, it is worth pointing out again that since shocks are not the only weak-solution discontinuities that arise in multimaterial MHD simulations, they are not the only locations where potential convergence failures are corrected by the new algorithm. In particular, contact discontinuities, material discontinuities, and ionization fronts all present potential problems for the naive Biermann algorithm, since they all represent locations where Pe changes discontinuously (despite the continuity of total pressure P at such locations), and are therefore all sites where the Biermann effect can be expected to generate a magnetic field. They all represent physical situations in which the new algorithm provides a correct treatment of the Biermann effect.

We would like to thank Gianluca Gregori for discussions that helped to focus this work on the relevant physical issues. We would also like to thank the anonymous referee for encouraging us to discuss more fully the astrophysical implications of this work and Andrey Kravtsov, Brian O'Shea, and John ZuHone for correspondence and conversations about the application of this work to the creation of magnetic fields by the Biermann effect in galaxy clusters. This work was supported in part by the National Science Foundation under grant AST-0909132, and in part by the U.S. DOE NNSA ASC through the Argonne Institute for Computing in Science under field work proposal 57789. The software used in this work was developed in part by funding from the U.S. DOE NNSA-ASC and OS-OASCR to the Flash Center for Computational Science at the University of Chicago.

APPENDIX A: SHOCK DETECTION

We describe here our shock detection algorithm, which we used in the verification of the Biermann effect at shocks. This algorithm has some benefits over shock-detection algorithms such as the ones described in Balsara & Spicer (1999), Miniati et al. (2000), and Ryu et al. (2003), in that it furnishes an estimate of the shock speed directly from a single time slice of spatial data; weights the mass, momentum, and energy components of the R–H relations equally; yields a single-cell-wide shock interface; verifies that local characteristics converge on the shock surface; and does not impose arbitrary thresholds on physical quantities such as compression ratios or Mach numbers.

We start with a mesh of cells containing fluid values (for present purposes it is immaterial whether these are cell averages or point values). The algorithm seeks a set of "shocked" cells satisfying the following criteria:

  • 1.  
    the shocked surface is one cell wide: each shocked cell is the location of a maximum of $|\nabla P|$ in the direction of $\nabla P$, where P is the total fluid pressure;
  • 2.  
    fluid quantities change correctly: the neighborhood of each shocked cell exhibits compression, acceleration, and heating in the direction of $\nabla P$;
  • 3.  
    R–H conditions obtain: the full R–H conditions on mass, momentum, and energy hold between the fluid upstream and downstream of each shocked cell, where "upstream" means in the direction ${\boldsymbol{n}} \equiv -\nabla P/|\nabla P|$;
  • 4.  
    characteristics converge on the shock surface: the characteristic convergence criterion holds at each shocked cell—the shock is supersonic upstream and subsonic downstream.

The result is an efficient and reliable algorithm that, among other things, yields an accurate estimate of the shock speed, which is essential to our verification work on the Biermann effect. Note that we neglect the induction equation in the shock detection algorithm and operate using only the material pressure—not the total pressure, which includes the magnetic pressure. For weakly magnetized plasmas, such as the ones we consider in this paper, this creates no difficulty in identifying the correct shock surface. Some generalization would be required for significantly magnetized plasmas, particularly to the third and fourth conditions above.

We now describe in somewhat greater detail the algorithm outlined above.

The Shocked Surface is One Cell Wide

The first condition above is tantamount to insisting that the shocked cell be sited at a position where the gradient of P is steepest. In general, volume-based hydrodynamic advance schemes spread shocks over several cells. If we were not to impose this condition, we would obtain a thickened shell of shocked cells, which would complicate the criterion (specified below) for classifying neighboring cells as adjoining the shock upstream or downstream.

The condition is very easy to enforce: in the immediate neighborhood of the cell being tested for "shockedness," we check the adjoining cells closest to the direction of the normal vector ${\boldsymbol{n}} $ and its mirror image $-{\boldsymbol{n}} $. This is done by forming the vector of integers ${\Delta }{\boldsymbol{i}} \equiv {\rm NINT}[{\boldsymbol{n}} /{\rm max} ({{n}_{1}},{{n}_{2}},{{n}_{3}})]$ (where the NINT function ascribes the nearest integer to each component of its real vector argument) and adding it to the cell discrete index vector ${\boldsymbol{i}} $, to reach the cells at ${\boldsymbol{i}} \pm {\Delta }{\boldsymbol{i}} $. If a stencil-based estimate of $|\nabla P|$ is greater in the candidate cell than in the two so-chosen adjoining cells, then the condition is passed successfully.

Fluid Quantities Change Correctly

An easy sanity check for "shockedness" of a cell is that cells downstream should unfailingly be (a) compressed, (b) accelerated (in the $+{\boldsymbol{n}} $-direction), and (c) heated, with respect to cells upstream.

We introduce a shock-width parameter h, such that the shock-adjoining cells "upstream" and "downstream" of the shocked cell are respectively located at $\pm h\times {\Delta }\times {\boldsymbol{n}} $ relative to the cell under study, where Δ is the length of a cell side. For FLASH with the HLL Riemann solver (Toro 2009) and second-order reconstruction, an appropriate value of h is 1.5. We then simply verify that the fluid variables ρ, ${\boldsymbol{u}} $, and eT (the total thermal energy) satisfy the conditions $\left[ \rho \right]_{u}^{d}\gt 0$, $\left[ {\boldsymbol{u}} \cdot {\boldsymbol{n}} \right]_{u}^{d}\gt 0$, and $\left[ {{e}_{T}} \right]_{u}^{d}\gt 0$.

R–H Conditions Obtain

The R–H conditions, which express flux conservation in the frame of the shock, have the form

Equation (A1)

Here ${\bf \Phi }$ is the state vector of conserved field quantities, ${{{\bf \Phi }}^{T}}\equiv [\rho ,\rho {\boldsymbol{u}} \cdot {\boldsymbol{n}} ,\rho \mathcal{E}]$, that is, the mass, normal momentum, and energy densities. ${\boldsymbol{F}} $ denotes the vector of fluxes corresponding to ${\bf \Phi }$. The subscripts "d" and "u" denote "downstream" and "upstream" states, respectively (Toro 2009).

We fit R–H conditions to the upstream/downstream data defined using the shock-width parameter h above. In this fit, the shock speed D is a free parameter to be adjusted to minimize a fit residual. Given a positive-definite inner product $\left[ {{{\Phi }}_{1}},{{{\Phi }}_{2}} \right]$ on the vector state space in Equation (A1), we may define the normalized residual R2 as

Equation (A2)

where ${\Delta }{\bf \Phi }\equiv {{{\bf \Phi }}_{d}}-{{{\bf \Phi }}_{u}}$ and ${\Delta }{\boldsymbol{F}} \equiv {\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm d}}})-{\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm u}}})$. We may then easily minimize R2 with respect to D, obtaining

Equation (A3)

Equation (A4)

Given these definitions, we consider that the R–H conditions are satisfied if R2Min is less than some small threshold value.

It remains to define the inner product $[\cdot ,\cdot ]$ used in these expressions. It is obvious that the naive unweighted sum-of-squares of the components of vectors such as ${\Delta }{\bf \Phi }$ is not an acceptable inner product, as it is dimensionally senseless. We require at a minimum a positive-definite metric that brings all vector components to the same physical dimensions, so that we do not wind up adding mass densities to energy densities, and so on. An additional desirable requirement is that all three components of Equation (A1) should contribute similar magnitudes to Equations (A3) and (A4), so that they are all weighted equally in the fit.

A simple metric that accomplishes these requirements may be built using the local sound speed cs of the candidate shocked cell:

Equation (A5)

where

Equation (A6)

This choice produces compatible dimensions and comparable magnitudes, because in and downstream of the shock we have ${{c}_{s}}\sim {{({{k}_{B}}T/{{m}_{i}})}^{1/2}}\sim D$, so that in the frame of the shock, the mass flux is $\rho D\sim \rho {{c}_{s}}$, the momentum flux has magnitude $P\sim \rho c_{s}^{2}$, and the energy flux has magnitude $\rho {{D}^{3}}\sim \rho c_{s}^{3}$.

We adopt this metric in computing Equations (A3) and (A4). We find that a threshold of $R_{{\rm Min}}^{2}\lt 0.01$ is adequate for identifying cells satisfying the R–H conditions with few false positives and false negatives, in the simulations exhibited in this paper.

Characteristics Converge on the Shock Surface

Our final criterion is easily stated and checked: shocks are supersonic upstream and subsonic downstream, so that characteristics from the family implicated in the shock converge on the shock (see Courant & Friedrichs 1948). This is an essential stability criterion for the solution. Using the shock speed DMin calculated while fitting the R–H conditions to the data, this condition is simply

Equation (A7)

where ${{c}_{s,d}}$ and ${{c}_{s,u}}$ are the downstream and upstream sound speeds, respectively, as determined with respect to the shock width h defined above.

Footnotes

Please wait… references are loading.
10.1088/0004-637X/802/1/43